NVU dynamics. II. Comparing to four other dynamics 
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In the companion paper [Ingebrigtsen et al.] an algorithm was developed for tracing out a geodesic 
curve on the constant-potential-energy hypersurface. Here simulations of this NVU dynamics are 
compared to results for four other dynamics, both deterministic and stochastic. First, NVU dy- 
namics is compared to the standard energy-conserving Newtonian NVE dynamics by simulations 
of the Kob-Andersen binary Lennard- Jones liquid, its WCA version (i.e., with cut-off's at the 
pair potential minima), and the Gaussian Lennard- Jones liquid. We find identical results for all 
quantities probed: radial distribution functions, incoherent intermediate scattering functions, and 
mean-square displacement as function of time. Arguments are then presented for the equivalence of 
NVU and NVE dynamics in the thermodynamic limit; in particular to leading order in 1/N these 
two dynamics give identical time-autocorrelation functions. In the final section NVU dynamics is 
compared to Monte Carlo dynamics, to a diffusive dynamics of small-step random walks on the 
constant-potential-energy hypersurface, and to Nose-Hoover NVT dynamics. If time is scaled for 
the two stochastic dynamics to make their single-particle diffusion constants identical to those of 
NVE dynamics, the simulations show that all five dynamics are equivalent at low temperatures 
except at short times. 



I. INTRODUCTION 

In the companion paper (Paper I [lj) we developed 
a stable numerical algorithm for tracing out a geodesic 
curve on the constant-potential-energy hypersurface f2 
of a system of N classical particles. If U(ri, rjy) is 
the potential energy as a function of the particle coordi- 
nates, for a given value Uq of the potential energy Q is 
the 3N — 1 dimensional Riemannian diffcrentiable mani- 
fold defined by (where R = (ri, ...,rjv) is the position in 
the 3JV dimensional configuration space) 



{R e R iN | U(R) = U } . 



(1) 



Geodesic motion on f2 is termed NVU dynamics in anal- 
ogy with standard Newtonian NVE dynamics, which 
conserves the total energy E. Motivations for studying 
NVU dynamics were given in Paper I. The present paper 
compares NVU dynamics to four other dynamics, two 
deterministic and two stochastic, concluding that NVU 
dynamics is a fully valid molecular dynamics. 

The path of shortest distance between two points on 
a Riemannian manifold is a so-called geodesic curve. By 
definition a geodesic is as a curve of stationary length, 
i.e., one for which small curve variations keeping the 
two end points R^ and Rb fixed to lowest order do not 
change the curve length, i.e., 



5 dl = 0. 



(2) 



By discretizing this condition and carrying out the vari- 
ation, keeping the potential energy fixed by introducing 
Lagrangian multipliers, a "basic NVU algorithm" was 
derived in Paper I (F is the 3A^-dimensional force vector 
and i is the time-step index): 



R?; 
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Fi • (Ri — Rj_i) 



(3) 



This algorithm works well, but for very long simulations 
numerical errors accumulate and U drifts to higher values 
("entropic drift", see Paper I). This problem is also expe- 
rienced for the total energy in NVE algorithms [2|, and 
it is not more severe for NVU than for NVE dynamics. 
A fully stable NVU algorithm was developed in Paper 
I, which may be summarized as follows. If one switches 
to the leap-frog representation and defines the position 
changes by A i+1 / 2 = Rj+i — Ri, the stable NVU al- 



gorithm is: A i+1 / 2 = 'oA,_ 1/2/ 
step length and A^_i/ 2 = 



1-1/2 



where Iq is the 



i-1/2 



-2F,- • A,- 
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Ui-i — J7o)Fj/Ff. Just as for standard NVE dynam- 
ics a final stabilization introduced is to adjust the po- 
sition changes slightly, e.g., every 100th step, in order 
to eliminate numerical drift of the center of mass coor- 
dinate. In the simulations reported below we used the 
stabilized NVU algorithm. However, since the stabiliza- 
tion is merely a technicality, the basic NVU algorithm 
Eq. ([3]) is used for theoretical considerations. 

Constant-potential-energy algorithms were previously 
considered in papers dating back to 1986 by Cotterill and 
Madsen et al. [3] and in 2002 by Scala et al. 4]. In the 
same spirit, but a slightly different context, Stratt and 
coworkers in 2007 and 2010 considered geodesic motion 
in the space of points with potential energy lower than Uq 
||. In the thermodynamic limit these points are almost 
all of energy very close to Uq. We refer to Paper I for a 
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more detailed discussion of how NVU dynamics relates 
to these earlier works. 

NVU dynamics invites to an alternative view of 
molecular motion. Instead of focusing on the stan- 
dard potential-energy landscape in (37V + 1) dimensions 
@, NVU dynamics adopts the configuration-space mi- 
crocanonical viewpoint and focuses on the (3N — 1)- 
dimensional Riemannian hypersurface f2. The classical 
potential-energy landscape picture draws attention to the 
stationary points of the potential-energy function, in par- 
ticular its minima, the so-called inherent states [6(. In 
contrast, all points on il have the same probability in 
NVU dynamics and there are no energy barriers - all 
barriers are of entropic nature defining unlikely parts of 
il that must be passed Despite the absence of 

energy barriers in the ordinary sense of this term, NVU 
dynamics is fully able to describe locally activated events 
(hopping processes between local potential-energy min- 
ima). The NVU "heat bath" is provided by the multi- 
tude of configurational degrees of freedom [3|-|5[ . 

The present paper compares NVU dynamics to other 
molecular dynamics, including stochastic ones. We first 
compare to NVE dynamics, which is also determinis- 
tic, and conclude that for large systems the two dy- 
namics are equivalent. We proceed to compare to other 
kinds of dynamics, inspired by previous works: The first 
investigation providing long-time simulations that com- 
pared different dynamics (Newtonian versus Langevin) 
was presented by Gleim et al. 7j. They studied the 
Kob-Andersen binary Lennard- Jones (KABLJ) mixture 
Q at different temperatures and found that below a cer- 
tain temperature (T < 0.8), the temperature dependence 
of the diffusion constant and of the structural relaxation 
time was identical in the two dynamics. This type of in- 
vestigation was extended by Szamel et al. [i| to Brownian 
dynamics, i.e., stochastic dynamics without the momen- 
tum degrees of freedom. They found power-law fitting 
exponents for the temperature dependence of the diffu- 
sion constant and relaxation time very close to those of 
NVE dynamics. Subsequently Berthier et al. in- 
vestigated Monte Carlo dynamics for which agreement 
with Newtonian dynamics was also established, both for 
a strong and a fragile model glass former (an Si02 model 
and the KABLJ model). This, however, did not apply 
for higher-order time-correlation functions, a fact con- 
tributed to the presence of different conservation laws 

We compare below NVU dynamics to four other 
dynamics: Newtonian dynamics (NVE), Nose- Hoover 
NVT dynamics (n), Monte Carlo dynamics (MC) [i~2j . 
and a diffusive small-step random-walk dynamics on the 
potential energy hypersurface (RW). Section |TT] com- 
pares NVU dynamics with the "true" (NVE) time evo- 
lution defined by Newton's second law. This is done 
by simulations of the KABLJ liquid, as well as of the 
Weeks-Chandler- Andersen (WCA) approximation [13j to 
the KABLJ liquid (KABWCA). SectionQIDgives some in- 
tuitive arguments for the equivalence of NVU and NVE 



dynamics in the thermodynamic limit. Section [IV] com- 
pares NVU dynamics with NVT, MC, and RW dynam- 
ics. Section fVl gives a brief summary and outlook. 



II. SIMULATIONS COMPARING NVU 
DYNAMICS TO NVE DYNAMICS 
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FIG. 1: (a) Probability density of (AU^nvu) 2 given by Eq. 
{5P for the Kob-Andersen binary Lennard-Jones (KABLJ) 
mixture at p = 1.2 and T = 0.44 (b) Probability density 
for (Ati !N vu) 2 - ((AU^Nvu) 2 } for 256, 1024, and 8192 parti- 
cles of the single-component LJ liquid (T — 0.70, p = 0.85), 
showing a narrowing as the particle number increases. 

In NVU dynamics a geodesic is traced in configuration 
space. Physically, this curve may be traversed with any 
velocity; comparing however to NVE dynamics suggests 
an obvious time measure for NVU dynamics, as we shall 
see now. Limiting ourselves for simplicity to systems of 
particles with identical masses m, the Verlet algorithm 
for NVE dynamics with time step At^vE is @,II3| 



R-i-l 



2 Ri — R,_i 



(A(we) 2 , 



(4) 



Comparing to Eq. ([3]) suggests the following identifica- 
tion of an NVU time step AUnvu 



(Ati^vu) 2 



Fj • (Rj — Rj-i) 



(5) 
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This quantity is identical to Iq\ of Paper I; our simula- 
tions show that the right-hand side is always positive for 
small Iq- We have no proof of this, but presumably it 
applies rigorously in the thermodynamic limit. 

In the following data are given in terms of the nat- 
ural units for the Lennard- Jones pair potential; for the 
KABLJ and KAWCA system length and energy are given 
in units of the large-particle parameters o~aa and caa, re- 
spectively. The simulation sizes are N = 1024, 1000, and 
1024 for KABLJ, KAWCA, and Lennard- Jones Gaussian, 
respectively (see below). 

The probability distribution of (A^jvvc/) 2 is given in 
Fig. [TJ for an N = 1024 KABLJ liquid at p = 1.2 and 
T = 0.44 [l5l |. The simulations behind this, as well as all 
below figures, were initiated by choosing the two initial 
configurations from a well-equilibrated NVE simulation. 
The potential energy U in the NVU simulation was cho- 
sen as U = (U)nve at the relevant state points. The 
probability distribution of Fig. [TJ is a Gaussian, which is 
consistent with the fact that (At^jyvi/) is a sum of many 

dons. 





FIG. 2: The radial distribution functions for the KABLJ sys- 
tem at p = 1.2. The black lines give results from NVE sim- 
ulations, colored circles from NVU simulation where green, 
red, and blue denote, respectively, AB, AA, and BB pairs, 
(a) T = 2.0; (b) T — 0.405. 



In view of the above, for comparing NVU and NVE 
generated sequences we define the NVU time step length 
Ai/vvt/ as the average of Eq. ([S]), i.e., 



(AtNvu) 2 



= -2 



Fj • (Rj — R-i-i) 



(6) 



First, we present results that compare static averages 
of NVU and NVE simulations. Figure [5] shows the three 
radial distribution functions for the KABLJ liquid at two 
different state points. Clearly the two algorithms give 
identical results. Figure [3] shows NVU and NVE results 
for the mean-square displacement and the incoherent in- 
termediate scattering function of the KABLJ liquid at 
density p = 1.2 over a range of temperatures. The mean- 
square displacement and the incoherent scattering func- 





FIG. 3: (a) Mean-square displacement and (b) inco- 
herent intermediate scattering function at the wave vec- 
tor of the first peak of the AA structure factor. Both 
simulations were performed at p = 1.2 and T = 
2.0, 0.80, 0.60, 0.50, 0.44, 0.42, 0.405 (left to right) for the 
KABLJ mixture (1024 particles). NVE dynamics is given 
by the filled black circles connected by straight lines, NVU 
dynamics by the red crosses. 



Corresponding figures are shown in Fig. H] for the 
Weeks-Chandler- Andersen (WCA) approximation, which 
cuts off interactions beyond the energy minima, i.e., keep 
only the repulsive part of the potential. The WCA 
version of the system has a similar structure, but a 
much faster dynamics in the supercooled regime [l6l . [l7| . 
Again, NVU and NVE dynamics give identical results. 
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change (A 2 A(t)) (i.e., the relative deviations go to zero 
as N — > oo). Consider the time-autocorrelation function 
of an extensive quantity A with zero average. In this case 
the time-autocorrelation function scales in both ensem- 
bles as N, and the proposed equivalence of the dynamics 
means that \(A(0)A(t)) NVU - (A(0)A(t)) NVE \ oc N° as 
N — > oo. Intuitively, what happens is that, since in NVE 
dynamics the relative potential-energy fluctuations go to 
zero N — > oo, it becomes a better and better approxima- 
tion 




FIG. 4: (a) Mean-square displacement and (b) incoherent 
intermediate scattering function at the same wave vector as 
in Fig. \3\ Both simulations were performed at p = 1.2 and 
T = 2.0, 0.80, 0.60, 0.50, 0.44 and 0.40 (left to right) for the 
WCA approximation to the KABLJ mixture. NVE dynamics 
is given by the hlled black circles connected by straight lines, 
NVU dynamics by the red crosses. 

We studied also the so-called Lennard- Jones Gaussian 
system defined by a pair potential that adds a Gaussian 
to a LJ potential [lH , a liquid that is not strongly corre- 



lating [16|. Figure [5] shows that also for this model the 



incoherent intermediate scattering function is the same 
for NVU and NVE dynamics. In conclusion, for all sys- 
tems simulated we found NVU = NVE. This applies 
even for N = 65 particles of the KABLJ liquid (T = 0.8, 
p=1.2). 
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FIG. 5: The incoherent intermediate scattering function at 
p = 0.8 and T — 1.4 for a Lennard-Jones Gaussian system 
\13] . The black circles represent an NVE simulation, the red 
symbols an NVU simulation. It should be noted that the 
system phase separated during the simulation. 



There exists in analytical mechanics a variational prin- 
ciple that does not involve time at all. This is the Mau- 
pertuis principle from 1746 [20L l2lj . a variational prin- 
ciple that is originally due to Jacobi and for this reason 
is sometimes referred to as "Jacobi's form of the least 
action principle" [l9|, [2lj|. This states that a classical- 
mechanical system of fixed energy E follows a curve in 
configuration space obeying (with fixed end points) 



/ v / 2m(E-U) dl = 0. 
Jr a 



(7) 



III. INTUITIVE ARGUMENTS FOR THE 
EQUIVALENCE OF NVU AND NVE DYNAMICS 

AS N ->• oo 

The above results raise the question: Are NVU and 
NVE dynamics mathematically equivalent in some well- 
defined sense? The two algorithms are not identical, of 
course; that would require no variation in the quantity 
&ti,NVU (Fig. [J). On the other hand, the Ati^vu dis- 
tribution narrows as the particle number increases (Fig. 
[T](b)). From this NVU and NVE dynamics are expected 
to become equivalent for N — > oo in the following sense: 
For any configurational quantity A, to leading order in 
1 /N there is identity of dynamic quantities like the time- 
autocorrelation function (A(0)A(i)) or the mean-square 



One may argue that the relative variations of the inte- 
grand go to zero as N — » oo. Thus the integrand in 
this limit becomes effectively constant and can be taken 
outside the variation, implying Eq. ([2]) for motion that 
effectively takes place on the constant potential energy 
surface [5j. 

If I is the path length parametrizing the path, Eq. ([7]) 
implies 0, 13 d 2 R/dl 2 = [F - (F ■ t)t]/2(£ - C/(R)) 
where t = dTL/dl is the unit vector tangential to the 
path. The term F — (F • t)t is the (vector) component 
of the force normal to the path. In the thermodynamic 
limit the path approaches the constant-potential-energy 
hypersurface O, i.e., F • t = 0. In this limit one has dl oc 
dt because the relative kinetic energy fluctuations go to 
zero. In this way the Maupertuis principle is equivalent 
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to both the geodesic equation Eq. ([2]) and to Newton's 
second law R = F/m. 
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FIG. 6: (a) The dynamical fluctuations quantified by xt(t) 
for the A particles at p — 1.2 for a system with 1024 parti- 
cles. The black circles give results for an NVE simulation, 
the red, green, and blue symbols represent NVU simulations 
at, respectively, T = 0.44, 0.42, 0.405. (b) The dynamical 
fluctuations quantified by %^{t) f° r ^ ne A particles at p — 1.2 
for a system with 2048 particles. The black circles give results 
for an NVE simulation, the violet and red symbols represent 
NVU simulations of, respectively, T = 0.50, 0.44, 0.42. In- 
creasing the number of particles does not appear to decrease 
the deviation between the two dynamics. 



The equivalence of NVU and NVE dynamics in the 
thermodynamic limit relates to static averages as well as 
to time-autocorrelation functions of extensive quantities. 
Just as one must be careful when comparing fluctuations 
between different ensembles, fluctuations relating to the 
dynamics need not be the same for NVU and NVE dy- 
namics. As an example, Fig. [6] shows the quantity X4,(i) 
defined by X4 (t) = N A [ (F 2 A (k, t)) - (F sA (k, t)) 2 ] for the 
KABLJ system at three temperatures and two values of 
N. Xi quantifies the incoherent intermediate scattering 
function fluctuations (Hf. For Xi(t) NVU and NVE dy- 
namics do not appear to give identical results. A related 
observation was made by Berthier et al. , who showed that 
X4 (t) is not the same in NVE and in NVT dynamics ■ 



IV. COMPARING NVU DYNAMICS TO NVT, 
MONTE CARLO, AND DIFFUSION ON Q 

This section compares simulations using NVU dynam- 
ics to results for three further dynamics, two of which are 
standard. We focus on the viscous regime. One dynamics 
is the Nose-Hoover NVT dynamics, a deterministic sam- 
pling of the NVT canonical ensemble that may be de- 
rived from a "virtual" Hamiltonian [n], [23| . The second 
standard dynamics considered is the Metropolis Monte 
Carlo (MC) algorithm, which generates a stochastic se- 
quence of states giving the correct NVT canonical en- 
semble distribution. The third dynamics employed below 
is also stochastic, it simulates diffusion on the constant- 
potential-energy hypersurface f2 by a small step-length 
random walk (RW) on SI. This was discussed by Scala et 
al. [3] , who proposed the following equation of motion 



where Ar/i is a 37V-dimensional random vector (see be- 
low). Equation © implies Fj R^ =0, which ensures the 
potential-energy conservation required for staying on fi. 

The RW algorithm was discretized and implemented 
as a " predictor-corrector" algorithm in the following way. 
A vector Arji was chosen from a cube with length L = 
O.Oler. This is small enough to ensure that the dynamics 
generates the correct NVE radial distribution function 
and at the same time has no effect on the average dy- 
namical quantities. Positions were updated via 



AtAn • F 

Ri+i = R t + AtA Vl -g ^F, . (9) 

F i 

Finally, Rj+i was corrected by applying two iterations of 
Rj+i = R,:+i — u% + 1 i ^ Ua Fi+i in order to eliminate long- 
time entropic drift of the potential energy. 

MC and RW dynamics involve no generic measures 
of time. We compared their results to NVU dynamics 
by proceeding as follows. At any given state point the 
time-scaling factor was determined from the long-time 
behavior of the mean-square displacement by requiring 
that the single-particle displacement obeys (Ax 2 (t)} = 
2Dt for t — > oo with the NVE diffusion constant D. By 
construction this ensures agreement with the long-time 
mean-square displacement of NVE dynamics. 

In Fig. [7] we show the incoherent intermediate scat- 
tering function of the KABLJ liquid for all investigated 
dynamics at several state points. 
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FIG. 7: The incoherent intermediate scattering function for 
all five investigated dynamics for the KABLJ mixture at p = 
1.2 and T = 2.0, 0.80, 0.60, 0.50 and 0.44. The black curve 
(with hlled circles) is the NVE simulation, red crosses: NVU , 
green squares: NVT, magenta diamonds: MC, blue triangles: 
RW. 



V. SUMMARY AND OUTLOOK 



NVU dynamics traces out geodesic curves on the 
3iV — 1 dimensional potential-energy hypersurface f2. 
We have compared NVU dynamics with other dynam- 
ics. Simulations supplemented by non-rigorous analyt- 
ical arguments showed that NVU and NVE dynamics 
are equivalent in the thermodynamic limit, i.e., typical 
autocorrelation functions become identical to leading or- 
der in 1/N. Furthermore, NVU dynamics was compared 
to two stochastic dynamics, standard Monte Carlo dy- 
namics and a small-step random walk on the constant- 
potential-energy hypersurface Q representing diffusion on 
n. Agreement was established for all dynamics, including 
also NVT dynamics, in the a-relaxation regime where 
inertial effects become unimportant. We conclude that 
NVU dynamics is a fully valid molecular dynamics. 



NVU and NVT dynamics agree quantitatively for all 
investigated state points. This is not surprising given 
the results of Sees. [TT] and Mil and the well-known 
fact that NVE and NVT dynamics give same time- 
autocorrelation functions to leading order in 1/N [2~4| . 
The incoherent intermediate scattering functions of MC 
and RW agree at all investigated temperatures. This is 
consistent with the recent results of Berthier et al. [l(| , 
who compared Langevin to MC dynamics. For lower 
temperatures (T < 0.80) quantitative agreement is found 
among all five dynamics investigated in the a-relaxation 
regime. 

The corresponding figure for the KABWCA system is 
shown in Fig. [5] The same conclusion is reached as for 
the KABLJ mixture. 
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FIG. 8: The incoherent intermediate scattering function for 
all five investigated dynamics for the KABWCA system at p = 
1.2 and T = 2.0, 0.80, 0.60, 0.50, 0.44 and 0.40. The black 
curve (with hlled circles) is the NVE simulation, red crosses: 
NVU, green squares: NVT, magenta diamonds: MC, blue 
triangles: RW. 



It is interesting to note that NVU dynamics, like 
any geodesic motion on a Riemannian manifold, can 
be formulated as a Hamiltonian dynamics based on the 
curved-space purely kinetic energy Hamiltonian H = 
1/2 J2 a b9 ab ( x )P a P h wnere x is the manifold coordinate, 
g a b the corresponding metric tensor, and p a the general- 
ized momenta (25[. Indeed, long ago Hertz argued that 
one should focus exclusively on the kinetic energy and 
describe classical mechanics as a geodesic motion on a 
high-dimensional Riemannian manifold (along the "ger- 
adeste Bahn" of this manifold, the straightest curve) [26| . 
Hertz' idea was to eliminate the force and potential en- 
ergy concepts entirely from mechanics and replace parti- 
cle interactions by constraints among the coordinates; the 
relevant manifold is defined by these constraints. This 
is not what we have done here. There is, however, the 
fundamental similarity between the Hertz and the NVU 
approaches that both build on the conceptual simplifi- 
cation of "replacing Newton's second law by Newton's 
first law". Moreover, as shown in Appendix A the ef- 
fect of masses enters into the metric of the Riemannian 
manifold in precisely the same way as we need for NVU 
dynamics when it is generalized to deal with systems of 
varying masses. Thus NVU dynamics realizes Hertz's 
ideas to a large extent. 

From a technical point of view NVU dynamics pro- 
vides few advantages because it is not faster than NVE 
or NVT dynamics. However, by referring directly to the 
mathematical properties of a Riemannian differentiable 
manifold NVU dynamics leads to a new way of thinking 
about the classical mechanics of many-particle systems. 
Future work should focus on relating the mathematical 
properties of to the physical properties of the system 
in question. It is our hope that in this way new insights 
into liquid dynamics may be arrived at by adopting the 
NVU viewpoint. 
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Appendix A: Generalization of the NVU algorithm 
to deal with systems of different particles masses 

Papers I and II deal with systems of particles with 
identical mass m. The basic NVU algorithm Eq. ([3]), 
however, is well defined and works perfectly well for 
any general classical mechanical system. The algorithm 
traces out a geodesic on fi that is independent of the 
particles' masses, a geometrical path entirely determined 
from the function U(ri, r/v). Equation (|5j|, which en- 
sures NVU — NVE in the thermodynamic limit, only 
works if all particles have mass m. The question arises 
if a generalization of Eq. ([3]) is possible ensuring that 
NVU = NVE as N — > oo also for systems of particles 
with different masses. 

If the fc'th particle mass is m k , we seek to modify the 
basic NVU algorithm such that it for the k'th particle 
as JV — ¥ oo converges to (where r*- fc ' is the coordinate of 
the fc'th particle, F' fc ) the force on it, and subscript j is 
the time step index) 



„(fe) _ 9 _(fc) __(*) 



m k 



(Al) 



If the average mass is (m), we define reduced masses by 



m k 



(to) 



(A2) 



A geodesic is defined by giving the shortest distance be- 
tween any two of its close by points. In Paper I and in Eq. 



(J2J) of the present paper the distance measure is given by 
the standard Euclidian distance dl 2 — J2 k dr^ ■ dr^ k > . A 
change of the metric leads to different geodesies, 
sider the following metric: 



Con- 



dl 2 



k 



m k dr^ 



(A3) 



This is precisely the metric discussed by Hertz in his 
mechanics long ago (26j. In the "Hertzian" metric the 
discritezed path length used in deriving the NVU algo- 



rithm is (Paper I) Ylj 

is the time step index), 
becomes 



J2 k ™>k (rf^ - r^J (where j 
Thus the variational condition 



= 0. 



(A4) 



From this it follows via the ansatz of constant step length 
that 



„(*0 



= 2rf - rf\ 2\Fj • (R, - Rj-i)]Ff /(m^j 



,(*0 



" , 1 - - *j-i - ■ K"-j ~ *-\>-lfl*j I Vk^j) ■ 

(A5) 

This translates into Eq. (|A1[) for a suitably chosen 
At; likewise the relative fluctuations of the term 2[Fj • 
(R 3 — R 3 _i)]/F| go to zero in the thermodynamic limit 
(N -> oo), such that NVU = NVE in this limit. 

It is important to note that when systems of varying 
masses are considered, there is also the option of ignor- 
ing the masses as done by Stratt and coworkers [5| . This 
approach, which is consistent with, e.g., Brownian me- 
chanics, leads to a perfectly admissible NVU dynamics. 
However, it does not correspond to NVE dynamics in 
the same way as the above proposed "Hertzian" version 
of NVU dynamics with varying masses. 



[1] T. S. Ingebrigtsen et al., the previous paper. 

[2] M. P. Allen and D.J. Tildesley, Computer Simulation of 
Liquids (Oxford Science Publications, Oxford, 1987); D. 
Frenkel and B. Smit, Understanding Molecular Simula- 
tion (Academic, New York, 2002). 

[3] R. M. J. Cotterill, Phys. Rev. B 33, 262 (1986); R. M. 
J. Cotterill and J. U. Madsen, in Characterizing Complex 
Systems, Ed. H. Bohr (World Scientific, Singapore, 1990), 
p. 177; J. Li, E. Piatt, B. Waszkowycz, R. Cotterill, and 
B. Robson, Biophys. Chem. 43, 221 (1992); R. M. J. 
Cotterill and J. U. Madsen, J. Phys.: Condens. Matter 
18 , 6507 (2006). 

[4] A. Scala, L. Angelani, R. Di Leonardo, G. Ruocco, and 
F. Sciortino, Phil. Mag. B. 82, 151 (2002). 

[5] C. Wang and R. M. Stratt, J. Chem. Phys. 127, 224503 
(2007); ibid. 127, 224504 (2007); C. N. Nguyen and R. 



M. Stratt, J. Chem. Phys. 133, 124503 (2010). 
[6] M. Goldstein, J. Chem. Phys. 51, 3728 (1969); F. H. 

Stillinger, Science 267, 1935 (1995); F. Sciortino, J. Stat. 

Mech.: Theory Exp. 2005, 35 (2005); A. Heuer, J. Phys.: 

Condens. Matter 20, 373101 (2008). 
[7] T. Gleim, W. Kob, and K. Binder, Phys. Rev. Lett. 81, 

4404 (1998). 

[8] W. Kob and H. C. Andersen, Phys. Rev. E 51, 4626 

(1995); ibid. 52, 4134 (1995). 
[9] G. Szamel, and E. Flenner, Europhys. Lett. 67, 779 

(2004); E. Flenner, and G. Szamel, Phys. Rev. E. 72, 

011205 (2005). 

[10] L. Berthier, and W. Kob, J. Phys.: Condens. Matter 
19, 205130 (2007); L. Berthier, Phys. Rev. E 76, 011507 
(2007). 

[11] S. Nose, J. Chem. Phys. 81, 511 (1984); W. G. Hoover, 



8 



Phys. Rev. A 31, 1695 (1985). 

[12] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. 
H. Teller, and E. Teller, J. Chem. Phys. 21, 1087 (1953). 

[13] D. Weeks, D. Chandler, and H. C. Andersen, J. Chem. 
Phys. 54, 5237 (1971). 

[14] L. Verlet, Phys. Rev. 159, 98 (1967). 

[15] All simulations were performed using a molecular dynam- 
ics code optimized for NVIDIA graphics cards, which is 
available as open source code at http://rumd.org 

[16] U. R. Pedersen, N. P. Bailey, T. B. Schr0der, and J. C. 
Dyre, Phys. Rev. Lett. 100, 015701 (2008); U. R. Ped- 
ersen, T. Christensen, T. B. Schr0der, and J. C. Dyre, 
Phys. Rev. E 77, 011201 (2008); T. B. Schr0der, U. R. 
Pedersen, N. P. Bailey, S. Toxvaerd, and J. C. Dyre, 
Phys. Rev. E 80, 041502 (2009); N. P. Bailey, U. R. Ped- 
ersen, N. Gnan, T. B. Schr0der, and J. C. Dyre, J. Chem. 
Phys. 129, 184507 (2008); N. P. Bailey, U. R. Pedersen, 
N. Gnan, T. B. Schr0der, and J. C. Dyre, J. Chem. Phys. 
129, 184508 (2008); T. B. Schr0der, N. P. Bailey, U. R. 
Pedersen, N. Gnan, and J. C. Dyre, J. Chem. Phys. 131, 

234503 (2009); N. Gnan, T . B. Schr0der, U. R. Peder- 
sen, N. P. Bailey, and J. C. Dyre, J. Chem. Phys. 131, 

234504 (2009); N. Gnan, C. Maggi, T . B. Schr0der, and 
J. C. Dyre, Phys. Rev. Lett. 104, 125902 (2010); T . B. 
Schr0der, N. Gnan, U. R. Pedersen, N. P. Bailey, and J. 
C. Dyre, J. Chem. Phys. 134, 164505 (2011). 

[17] U. R. Pedersen, T. B. Schr0der, and J. C. Dyre, Phys. 
Rev. Lett. 105, 157801 (2010). 



[18] V. V. Hoang, and T. Odagaki, Physica B 403, 
3910 (2008). The Lennard- Jones Gaussian pair po- 
tential is e [(cr/r) 12 - 2(cr/r) 6 - e exp((r - r ) 2 /2ct 2 )] 
where ao = 0.14, eo = 1-5, ro = 1.47. 

[19] E. T. Whittaker, A treatise on the Analytical Dynamics 
of Particles and Rigid Bodies 4th Ed. (Cambridge Uni- 
versity Press, Cambridge, UK, 1999). H. Goldstein, Clas- 
sical Mechanics (Addison- Wesley, Reading, MA, 1950). 

[20] L. D. Landau and E. M. Lifshitz, Mechanics 2nd Ed. 
(Pergamon Press, Oxford, 1969). 

[21] Wikipedia article "Maupertuis' principle" 
( http: / /wikipedia. org ) . 

[22] C. Toninelli, M. Wyart, L. Berthier, G. Biroli, and J.-P. 
Bouchaud, Phys. Rev. E 71, 041505 (2005). 

[23] S. Toxvaerd, Mol. Phys. 72, 159 (1991); T. Ingebrigtsen, 
O. J. Heilmann, S. Toxvaerd, and J. C. Dyre, J. Chem. 
Phys. 132, 154106 (2010). 

[24] D. J. Evans, and B. L. Holian, J. Chem. Phys 83, 4069 
(1985). 

[25] Wikipedia article "Geodesies as Hamiltonian flows" 
( |http://w ikipedia.org[| . 

[26] H. Hertz, Die Prinzipien der Mechanik, in neuem Zusam- 
menhange dargestellt (Leipzig, 1894); J. Liitzen, Arch. 
Hist. Exact Sci. 49, 1 (1995); J. Liitzen, Mechanistic im- 
ages in geometric form: Heinrich Hertz's "Principles of 
mechanics" (Oxford University Press, Oxford, 2005); J. 
Preston, Stud. Hist. Phil. Sci. 39, 91 (2008). 



